package utils.math;

import java.util.ArrayList;
import java.util.List;

import misc.DerivativeFunction;
import misc.Function;
import misc.Seria;

public class MathCodes {

    /**
     * This method will calculate the value of given polynom at point x, using
     * Horner's algorithm.
     * 
     * @param polyCoefficient
     * @param x
     * @return The value of the polynom at point x
     */
    public static double horner(double[] polyCoefficient, double x) {
        return horner(polyCoefficient, x, 0);
    }

    /**
     * This method will calculate one panel, starting at x0, according to
     * Boole's Rule.
     * 
     * @param f
     * @param x0
     * @param h
     * @return
     */
    public static double calculatePanelAreaBR(Function f, double x0, double h) {
        double result = 0.0;

        result = (2 * h / 45)
                * (7 * f.valueAtPoint(x0) + 32 * f.valueAtPoint(x0 + h) + 12
                        * f.valueAtPoint(x0 + 2 * h) + 32
                        * f.valueAtPoint(x0 + 3 * h) + 7 * f.valueAtPoint(x0
                        + 4 * h));

        return result;
    }

    /**
     * This method will search a section (a,b) where f(x)=0
     * 
     * @param f
     * @param x0
     *            - first poin to start with.
     * @param h
     * @return an array of two double points - the section.
     */
    public static double[] blindSearch(Function f, double x0, double h) {
        double[] result = new double[2];
        double f1, f2;
        double a = x0, b = x0 + h;

        for (int k = 0; k < Configuration.MAX_TRIES; k++) {

            f1 = f.valueAtPoint(a);
            f2 = f.valueAtPoint(b);

            if (f1 * f2 < 0) {
                result[0] = Math.min(a, b);
                result[1] = Math.max(a, b);
                return result;
            }

            if (Math.abs(f1) < Math.abs(f2)) {
                a -= h * (a - b);
            } else {
                b += h * (b - a);
            }
        }
        result[0] = 0;
        result[1] = 0;

        return result;
    }

    /**
     * This method will calculate an Integral of function between (a,b) by using
     * Boole's Rule.
     * 
     * @param f
     * @param a
     * @param b
     * @param panelsCount
     * @return
     */
    public static double calculateIntegralBR(Function f, double a, double b, int panelsCount) {

        double result = 0.0;
        double panelSize = Math.abs(b - a) / (double) panelsCount;
        double h = panelSize / 4.0;
        for (int k = 0; k < panelsCount; k++) {
            result += calculatePanelAreaBR(f, a + (double)k * panelSize, h);
        }
        return result;
    }

    /**
     * Calculates Aitkein series.
     * 
     * 
     * @param seria
     * @param stopIndex
     * @param tolerance
     */
    public static List<List<Double>> calculateAitkeinSeries(Seria seria, int stopIndex, double tolerance ){
	
	List<Double> orgSeria = generateSeria(seria, stopIndex);
	List<Double> mistakes = new ArrayList<Double>();
	List<Double> accSeria = new ArrayList<Double>();
	
	for (int k = 2; k < stopIndex ; k++){
	    
	    Double mistake = -1 * (Math.pow(orgSeria.get(k-1)- orgSeria.get(k-2), 2)) /
		    		(orgSeria.get(k) - 2 * orgSeria.get(k-1) + orgSeria.get(k -2));
	    
	    mistakes.add(mistake);
	    
	    //calculate the accesselerated seria
	    accSeria.add(orgSeria.get(k-2) + mistake);
	    
	    if (Math.abs(mistake) < tolerance){
		
		break;
	    }
	}
	
	List<List<Double>> answer = new ArrayList<List<Double>>();
	answer.add(0, orgSeria);
	answer.add(1, mistakes);
	answer.add(2, accSeria);
	return answer;
    }
    
    
    /**
     * 
     * @param seria
     * @param stopIndex
     * @param s
     * @param tolerance
     * @return
     */
    public static List<List<Double>> calculateRichardsonSeries(Seria seria, int stopIndex,int s, double tolerance ){
	List<Double> orgSeria = generateSeria(seria, stopIndex);
	List<Double> mistakes = new ArrayList<Double>();
	List<Double> accSeria = new ArrayList<Double>();
	for (int k = 2; k <= stopIndex ; k = k + 2){
	    int j = k / 2;
	    Double mistake = (orgSeria.get(k) - orgSeria.get(j)) / (Math.pow(2, s) - 1);
	    mistakes.add(mistake);   	    
	    if (Math.abs(mistake) < tolerance){
		System.out.println("mistake < tolerance");
		break;
	    }
	}
	int counter = 0;
	for (Double eps : mistakes){
	    
	    accSeria.add(orgSeria.get(counter++ * 2) + eps );
	}
	List<List<Double>> answer = new ArrayList<List<Double>>();
	answer.add(0, orgSeria);
	answer.add(1, mistakes);
	answer.add(2, accSeria);
	return answer;
    }
    
    
    /**
     * 
     * @param f - function
     * @param a - initial section left limit 
     * @param b - initial section right limit
     * @param tolerance - cimputation tolerance
     * 
     */
    public static List<Double> bisection(Function f, double a, double b, double tolerance){
	
	double[] sSection = blindSearch(f, a, b);
	double x = sSection[0];
	double y = sSection[1];
	double middle = 0d;
	int row = 1;
	File.FileHandler.openFile("Section");
	File.FileHandler.openNewSheet("Section1");
	File.FileHandler.addNewLable("eps");
	File.FileHandler.addNewLable("(a + b) / 2");
	File.FileHandler.addNewLable("a");
	File.FileHandler.addNewLable("b");
	while( y - x > tolerance){
	    middle = (x + y) / 2.0;
	    double middleValueAtPoint = f.valueAtPoint(middle);
	    if (middleValueAtPoint * f.valueAtPoint(y) < 0) {
		x = middle;
	    } else if (f.valueAtPoint(x) * middleValueAtPoint < 0) {
		y = middle;
	    }
	    File.FileHandler.addNewFloatValue(0, row, y-x);
	    File.FileHandler.addNewFloatValue(1, row, middle);
	    File.FileHandler.addNewFloatValue(2, row, x);
	    File.FileHandler.addNewFloatValue(3, row, y);
	    row++;
	}
	File.FileHandler.closeFile();
	List<Double> answer = new ArrayList<Double>();
	answer.add(0, y - x);	 // mistake
	answer.add(1, middle);	// result
	return answer;
    }
    public static List<Double> bisection(final double[] polyCoef, double a, double b, double tolerance){
	
	Function f = new Function(){
	    @Override
	    public double valueAtPoint(double point) {
		return horner(polyCoef, point);
	    } 
	};
	return bisection(f,a,b,tolerance);
    }
    
    
    public static double NewtonRapson(DerivativeFunction f,double x, double y, double tolerance){
	
	double[] sSection = blindSearch(f, x, y);
	double a = sSection[0],  b = sSection[1], x0 = 0d;

	boolean val1 = ((f.getFirstDerivative().valueAtPoint(a) >= 0 & f.getFirstDerivative().valueAtPoint(b) >= 0)
			| (f.getFirstDerivative().valueAtPoint(a) < 0 & f.getFirstDerivative().valueAtPoint(b) < 0));	
	boolean val2 = ((f.getSecondDerivative().valueAtPoint(a) >= 0 & f.getSecondDerivative().valueAtPoint(b) >= 0)
			| (f.getSecondDerivative().valueAtPoint(a) < 0 & f.getSecondDerivative().valueAtPoint(b) < 0));
	if (val1 & val2){
	    if (f.getFirstDerivative().valueAtPoint(a) >= 0 & f.getSecondDerivative().valueAtPoint(a) <= 0){
		x0 = a;
	    }
	    if (f.getFirstDerivative().valueAtPoint(a) <= 0 & f.getSecondDerivative().valueAtPoint(a) >= 0){
		x0 = a;
	    }
	    if (f.getFirstDerivative().valueAtPoint(a) >= 0 & f.getSecondDerivative().valueAtPoint(a) >= 0){
		x0 = b;
	    }
	    if (f.getFirstDerivative().valueAtPoint(a) <= 0 & f.getSecondDerivative().valueAtPoint(a) <= 0){
		x0 = b;
	    }
	}else{	    
	    System.out.println("Derivative not Compatible in section ["+a+","+b+"]");
	    return 0d;
	}
	File.FileHandler.openFile("NR");
	File.FileHandler.openNewSheet("sheet 1");
	File.FileHandler.addNewLable("X");
	File.FileHandler.addNewLable("error");
	
	double x1 = 0d, error = 1d;
	int row = 2;
	File.FileHandler.addNewFloatValue(0, 1, x0);
	while (Math.abs(error) > tolerance){
	    x1 = x0 - (f.valueAtPoint(x0)) / f.getFirstDerivative().valueAtPoint(x0);
	    error = x1 - x0;
	    x0 = x1;
	    File.FileHandler.addNewFloatValue(0, row, x1);
	    File.FileHandler.addNewFloatValue(1, row, error);
	    row++;
	}
	File.FileHandler.closeFile();
	return x1 + error;
    }
    
    
    public static double secont(DerivativeFunction f,double x, double y, double tolerance){
	
	double[] sSection = blindSearch(f, x, y);
	double a = sSection[0], b = sSection[1], x0 = 0d, x1 = 0d;
	int row = 1;
	File.FileHandler.openFile("secont");
	File.FileHandler.openNewSheet("S1");
	File.FileHandler.addNewLable("X");
	File.FileHandler.addNewLable("error");
	boolean val1 = ((f.getFirstDerivative().valueAtPoint(a) >= 0 & f.getFirstDerivative().valueAtPoint(b) >= 0)
		| (f.getFirstDerivative().valueAtPoint(a) < 0 & f.getFirstDerivative().valueAtPoint(b) < 0));
        boolean val2 = ((f.getSecondDerivative().valueAtPoint(a) >= 0 & f.getSecondDerivative().valueAtPoint(b) >= 0)
        		| (f.getSecondDerivative().valueAtPoint(a) < 0 & f.getSecondDerivative().valueAtPoint(b) < 0));
        if (val1 & val2){
            if (f.getFirstDerivative().valueAtPoint(a) >= 0 & f.getSecondDerivative().valueAtPoint(a) <= 0){
        	x0 = a;
        	x1 = b;
            }
            if (f.getFirstDerivative().valueAtPoint(a) <= 0 & f.getSecondDerivative().valueAtPoint(a) >= 0){
        	x0 = a;
        	x1 = b;
            }
            if (f.getFirstDerivative().valueAtPoint(a) >= 0 & f.getSecondDerivative().valueAtPoint(a) >= 0){
        	x0 = b;
        	x1 = a;
            }
            if (f.getFirstDerivative().valueAtPoint(a) <= 0 & f.getSecondDerivative().valueAtPoint(a) <= 0){
        	x0 = b;
        	x1 = a;
            }
        }else{
            System.out.println("Derivative not Compatible in section ["+a+","+b+"]");
            return 0d;
        }
	
        double y0 = 0d, error = 1d;
	while (Math.abs(error) > tolerance){
	    y0 = x1 - ((x1 - x0) * f.valueAtPoint(x1) / ((f.valueAtPoint(x1) - f.valueAtPoint(x0) )));
	    error = y0 - x1;
	    x0 = x1;
	    x1 = y0;
	    File.FileHandler.addNewFloatValue(0, row, y0);
	    File.FileHandler.addNewFloatValue(1, row, Math.abs(error));
	    row++;
	}
	File.FileHandler.closeFile();
	return x1;
    }
    
    /**
     * Generates Seria elements. 
     * 
     * @param seria
     * @param stopIndex
     * @return
     */
    private static List<Double> generateSeria(Seria seria, int stopIndex) {
	
	List<Double> answer = new ArrayList<Double>();
	answer.add(0, new Double(0));
	for (int i = 1; i <= stopIndex; i++){
	    answer.add(i , seria.getElement(i));
	}
	
	return answer;
    }

    /**
     * 
     * @param polyCoefficient
     * @param x
     * @param index
     * @return
     */
    private static double horner(double[] polyCoefficient, double x, int index) {
        if (index == polyCoefficient.length - 1)
            return polyCoefficient[index];
        return polyCoefficient[index] + x
                * horner(polyCoefficient, x, index + 1);
    }
    
    
    

}